#include <eoprec_constdet_wilstype_fermact_w.h>
#include <actions/ferm/fermacts/eoprec_twm_fermact_array_w.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#define gaugeSiteSize 18 // real numbers per link
#include "chroma.h"
using namespace Chroma;
#include "tmf_4d_inv.h"
#include <complex>
#include "reunit.h"

void chroma_tmf_invert( LatticeFermion &eta, LatticeFermion &psi, multi1d<LatticeColorMatrix> &u,double M5);

/*
 * Test application that will make use of the 4d_dwf_invert class
 */
int main (int argc,char **argv )
{
    Chroma::initialize(&argc, &argv);

    /*
     * Read commandline arguments for test program
     */
    if ( argc != 5 )
    {
        printf("Usage: %s Lx Ly Lz Lt \n All must be even\n",argv[0]);
        exit(-1);
    }

    /*
     * Setting up the QDP enviroment
     */
    multi1d<int> nrow(Nd);
    nrow[0] = atoi(argv[1]);
    nrow[1] = atoi(argv[2]);
    nrow[2] = atoi(argv[3]);
    nrow[3] = atoi(argv[4]);

    Layout::setLattSize(nrow);
    Layout::create();

    /*
     * Create our '5th' dimension as array of fermions
     */
    LatticeFermion eta;
    LatticeFermion qudapsi;
    LatticeFermion chromapsi;
    
    eta.elem(all.start()).elem(0).elem(0).real() = 1;

    multi1d<LatticeColorMatrix> u(Nd);
    for ( int mu = 0; mu < Nd; mu++ )
    {
        //gaussian(u[mu]);
        //TEST_LIB::reunit(u[mu]);
        u[mu] = Real(1.0);
      //  HotSt(u);
    }

    initQuda(0);

    QDPIO::cout << "////////////////////////////////////////////\n";
    QDPIO::cout << "// Finding result using QUDA         \n";
    QDPIO::cout << "////////////////////////////////////////////\n";
    
    QDPIO::cout << "\n~Source Norm2 = " << norm2(eta) << " \n\n";

    QDPIO::cout << "Solving for the solution vector\n";
    tmf_4d_inv::twistedMassFermionInvert( eta, qudapsi, u, 1.8, 0);
    QDPIO::cout << "~Solution Norm2 = " << norm2(qudapsi) << " \n\n";


    QDPIO::cout << "////////////////////////////////////////////\n";
    QDPIO::cout << "// Finding result using CHROMA                 \n";
    QDPIO::cout << "////////////////////////////////////////////\n";

    QDPIO::cout << "\n~Source Norm2 = " << norm2(eta) << " \n\n";

    chroma_tmf_invert(eta, chromapsi, u, 1.8);
    QDPIO::cout << "~QDP Solution Norm2 = " << norm2(chromapsi) << " \n\n";

    QDPIO::cout << "Comparing Results\n";
    QDPIO::cout << "~4dSolution - QDPSolution Norm2 = " << norm2(qudapsi-chromapsi) << " \n";


    endQuda();

    Chroma::finalize();
}




void chroma_tmf_invert( LatticeFermion &eta, LatticeFermion &psi, multi1d<LatticeColorMatrix> &u,double m5){
    /********************************************************
     * QDP Linop
     ********************************************************
     */
    typedef LatticeFermion T;
    typedef multi1d<LatticeColorMatrix> U;
    
    Real mass(m5);
    Real mq = 0;// mass;
    #define Ncg 3
    Real residuals[Ncg] ; 
    int  max_iter[Ncg];
    residuals[0] = 1.0e-6;
    residuals[1] = 1.0e-7;
    residuals[2] = 1.0e-8;
    //  residuals[3] = 1.0e-8;

    for(int i=0;i<Ncg;i++) max_iter[i]  = 1000;


    multi1d<int> bcs(Nd);
    bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;

    Handle< FermBC<T,U,U> > fbc(new SimpleFermBC< T, U, U >(bcs));
    Handle<CreateFermState<T,U,U> > cfs( new CreateSimpleFermState<T,U,U>(fbc));

    EvenOddPrecTwmFermActArrayParams twmparam;
    twmparam.Mass = mass;

    EvenOddPrecTwmFermActArray  S_f(cfs,twmparam );
    Handle< FermState<T,U,U> > fs( S_f.createState(u) );

    Handle< EvenOddPrecLinearOperatorArray<T,U,U> > M(S_f.precLinOp(fs,mq));
    Handle< LinearOperatorArray<T> > HM(S_f.precLinOp(fs,mq));

    // Check the result
    SysSolverCGParams invParam;
    invParam.RsdCG        = residuals[1];
    invParam.RsdCGRestart = residuals[1];
    invParam.MaxCG        = max_iter[1];
    invParam.MaxCGRestart = max_iter[1];

    Handle<LinOpSystemSolverArray<T> > Hsolver = new LinOpSysSolverCGArray<T>(HM,invParam);
    Handle<PrecFermAct5DQprop<T,U,U> > solver = new PrecFermAct5DQprop<T,U,U>(M,Hsolver);
    SystemSolverResults_t  stats =  (*solver)(psi,eta);
    int n_count = stats.n_count;

}
